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Abstract 

Density functional theory has been used to study the adsorption of molecular 
H2 on a graphene layer. Different adsorption sites on top of atoms, bonds 
and the center of carbon hexagons have been considered and compared. We 
conclude that the most stable configuration of H2 is physisorbed above the 
center of an hexagon. Barriers for classical diffusion are, however, very small. 

I. INTRODUCTION 

The adsorption of hydrogen by different forms of carbon has been studied by different 
groups. 0~! Dillon et al were the first to study the storage of molecular hydrogen by 
assemblies of single wall carbon nanotubes (SWCNT) and porous activated carbon. They 
pointed out that the attractive potential of the walls of the pores makes it possible a high 
density storage. From temperature-programmed desorption experiments Dillon et al con- 
cluded that those forms of carbon are promising candidates for hydrogen storage, although 
the density of hydrogen is still low in order to meet the requirements of the DOE Agency 
for novel hydrogen storage systems. More recently Levesque et a/i, Ye et a/i, and Liu et al 
I also studied the adsorption of molecular hydrogen on SWCNT at different temperatures 
and pressures. Chambers et al i have reported obtaining an extraordinary storage capacity 
by some graphite nanofibers but Wang and Johnson @ have tried unsuccessfully to confirm 
the high storage capacity by graphite nanofibers (slit pores) and SWCNT. Hynek et al i 
investigated ten carbon sorbents but only one of them could augment the capacity of com- 
pressed hydrogen gas storage vessels. The improvement was marginal at 190 K and 300 K 
but nonexistent at 80 K. The storage capacity of carbon nanotubes and graphitic fibers has 
been enhanced by doping with lithium and other alkali elements!. The alkali atoms seem 
to have a catalytic effect in dissociating the H 2 molecule and promoting atomic adsorption. 
An advantage is that the doped systems can operate at moderate temperatures and ambient 
pressure. 
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Some of the authors cited above Hi'i have also performed computer simulations of the ad- 
sorption of molecular hydrogen inside, outside and in the interstices of an array of SWCNT 
and in idealized carbon slit pores using model pair potentials to describe the interactions. 
Wang and Johnson @ adopted the semiempirical pair potential of Silvera and Goldman 
for the H2 — H2 interaction and the H2 — C interaction was modelled by a potential derived 
by Crowell and Brown by averaging an empirical Lennard- Jones H2 — C potential over 
a graphite plane. In the simulations Wang and Johnson used a hybrid path integral-Monte 
Carlo method. Johnson i also studied the influence of electrical charging of the tubes. Stan 
and Cole performed calculations based on a sum of isotropic Lennard- Jones interactions 
between the molecule and the C atoms of the tube. They calculated the adsorption poten- 
tial of a hydrogen molecule, considered as a spherically symmetric entity, as a function of 
distance from the axis of a SWCNT, along radial lines upon the center of an hexagon of 
carbon atoms and upon a carbon atom respectively. Those simulations give useful insight to 
interpret the results of the experiments. However the description of the interaction between 
H 2 and the graphitic surfaces of the SWCNT or the slit pores in those works is too simple. 
Simplicity is a necessary requirement for massive simulations involving several hundred (or 
several thousand) H2 molecules and an assembly of SWCNT of realistic size, but one can 
expect more realistic results if the interaction potential is derived from an ab initio calcu- 
lation. The adsorption of "atomic" hydrogen on a planar graphene sheet, that is a planar 
layer exfoliated from graphite, has been studied previously tlB Bercu and Grecu III used a 
semiempirical molecular orbital LCAO treatment at the INDO (intermediate neglect of dif- 
ferential overlap) unrestrictedHartree-Fock level and Jeloaica and Sidis 111 used the Density 
Functional formalism (DFT)ll3. In both works the description of the graphene layers was 
simplified by modelling this layer by a finite cluster C24-H42, where the hydrogen atoms sat- 
urate the dangling bonds on the periphery of the planar cluster. But, as mentioned above, 
hydrogen is adsorbed in molecular form by graphitic surfaces (SWCNT and slit pores), so 
in this work we study the interaction of an H2 molecule with a planar graphene layer. Since 
the graphene layers interact weakly in bulk graphite, the interaction of H2 with a graphitic 
surface is a localized phenomenon restricted to the outermost plane. For this reason our cal- 
culations have relevance for understanding the adsorption of H 2 on the walls of slit pores in 
graphite, and also for the case of adsorption by SWCNT, since these differ from a graphene 
layer only in the curvature of the layer. 



II. THEORETICAL METHOD AND TESTS 



To calculate the interaction between H2 and a planar graphene layer we use the ab 
initio fhi96md code, developed by Scheffler et al 0. This code uses the DFT 111 to compute 
the electronic density and the total energy of the system,and we have chosen the local 
density approximation (LDA) for exchange and correlation. 113 Only the two electrons of the 
H 2 molecule and the four external electrons (2s 2 2p 2 ) of each carbon atom are explicitly 
included in the calculations, while the Is 2 core of carbon is replaced by a pseudopotential. 
For this purpose we use the nonlocal norm-conserving pseudopotentials of Hamann et al 
0'0. Nonlocality in the pseudopotential is restricted to I = 2, and we take as local part of 
the pseudopotential the s component. The code employs a supercell geometry and a basis 
of plane waves to expand the electronic wave functions.^ 
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First we have tested the method for pure graphite. By minimization of the total energy 
with respect to the interatomic distances we obtained an in-plane C — C bond length equal 
to 2.66 a.u. and a distance between planar graphitic layers of 6.275 a.u. The corresponding 
experimental values ill are 2.68 a.u. and 6.34 a.u. respectively. The small (1%) underes- 
timation of bond lengths is characteristic of the LDA. Next we have studied an isolated 
graphene layer. Since the computer code uses a periodic supercell method, the cell axis has 
to be large in the z-direction to avoid the interaction between graphene sheets in different 
cells. Table I gives the calculated energy of the graphene layer as a function of the length c 
of the unit cell in the z direction, or in other words, as a function of the distance between 
parallel graphene layers. Results given for c = 20, 25, 30 and 35 a.u. show that the energy 
is well converged for those layer separations and that for c = 20 a.u. the error in the energy 
per atom is only about 1 in 10 5 . A cutoff of 40 Ry was used in all the calculations. We 
also tested the method by calculating the energy of the H 2 molecule, that was placed at the 
center of a simple cubic supercell. The total energy obtained for a plane wave cut-off energy 
of 40 Ry and supercell lattice constants of 18 a.u. and 20 a.u. is the same, -2.247 Ry as well 
as the bondlength, 1.48 a.u. Notice that this bond length is small compared to the C — C 
bond length. Anticipating the geometry to be used in the study of the interaction between 
H 2 and graphene, another set of calculations were performed for the energy of H 2 by placing 
the molecule in the superlattice described above in the study of the graphene layer, but this 
time without graphene. Calculations for distances between the imaginary graphene planes 
ranging from 20 a.u. to 35 a.u. (the plane wave cut-off was again 40 Ry) gave energies for 
the H 2 molecule, identical to the energies obtained for the cubic superlattice geometry. 



III. INTERACTION BETWEEN H 2 AND THE GRAPHENE LAYER 

For the periodicity of the system we have selected a unit cell with eight carbon atoms 
and one hydrogen molecule (see Fig. 1). If we place a hydrogen molecule at any point of the 
cell, the distance from this molecule to other in the nearest cells is 9.224 a.u. This separation 
is large compared to the bond length of H 2 (1.480 a.u.), and we have verified that there is no 
interaction between two hydrogen molecules separated by that distance. The interaction of 
the H 2 molecule with the graphene sheet has been studied by performing static calculations 
for two orientations of the axis of the molecule: axis perpendicular to the graphene plane and 
axis parallel to that plane. Three possibles configurations, called A, B and C below, have 
been selected for the perpendicular approach of the molecule to the plane: (A) upon one 
carbon atom, (B) upon the center of a carbon-carbon bond, and (C) upon the center of an 
hexagon of carbon atoms. On the other hand, for the parallel approach the molecule is placed 
upon the center of an hexagon of carbon atoms with the molecular axis perpendicular to two 
parallel sides of the hexagon, and this is called configuration D. These four configurations 
are given in the bottom panel of Fig. 1. To obtain the interaction energy curve for each of 
those four cases, the distance between the hydrogen molecule and the graphene layer was 
varied while maintaining the relative configuration. In these calculations the bond length of 
the H 2 was held fixed at 1.48 a.u., the bondlength of the free molecule. This is expected to 
be valid in the relevant region of the interaction. This constraint will, however, be relaxed 
in simulations described at the end of this section. Calculations were first performed in 
the parallel configuration (D) for a superlattice such that the distance between graphene 
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layers is 30 a.u. The plane wave cut-off was 40 Ry. The interaction energy curve is plotted 
in Figure 2 and the curve has a minimum at 5.07 a.u. For separations larger than this 
value the energy rises fast and reaches its asymptotic value for 10 - 11 a.u. The energy at 
the maximum possible separation between the center of mass of the H 2 molecule and the 
graphene plane for this superlattice, 15 a.u., was taken as the zero of energy. The figure also 
gives the results of a similar calculation for a smaller superlattice, such that the distance 
between graphene layers is 20 a.u. The corresponding energy curve, referred to the same zero 
of energy as above, is practically indistinguishable from the former curve. The calculations 
also show that for all practical purposes the energy curve has reached its asymptotic value 
for a distance of 10 a.u., that is the longest separation allowed for the superlattice of 20 a.u. 
This indicates that calculations using the smaller superlattice are enough for our purposes 
of studying the H 2 - graphene interaction. Then, the results of calculations corresponding to 
configurations A, B, C and D for a superlattice of 20 a.u. are given in Figure 3. The potential 
energy curves for the perpendicular approach (A, B, C) rapidly merge with each other for 
large H 2 - graphene separation, becoming indistinguishable from one another beyond 6.5 
a.u. Actually, curves A and B are very close in the whole range of separations although 
B is marginally more attractive. The common value of the energy of curves A, B and C 
at separation 10 a.u. is taken as zero of energy in Figure 3. The predicted equilibrium 
positions and the binding energies (depth of the minimum) of the different curves are given 
in Table II. The small magnitude of the binding energies, less than 0.1 eV, shows that 
the system is in the physisorption regime. Comparison of the four curves reveals that the 
most favorable position for the H 2 molecule is physisorbed in a position above the center 
of a carbon hexagon, and that the parallel configuration is slightly more favorable than the 
perpendicular one. We have verified that different orientations of the molecular axis with 
respect to the underlying carbon hexagon in the parallel configuration lead, in all cases, to 
the same curve D plotted in Fig. 3. The differences in binding energy shown in Table II are 
very small. For instance, configuration D and A only differ by 16 meV, and configuration D 
and C by 3 meV. 

Figures 4a and 4b give the electron density of the pure graphene layer in two parallel 
planes, 5 and 3 a.u. above the plane of the nuclei, respectively. The former one is very 
close to the preferred distance of approach for the H 2 molecule in configuration D. First of 
all one can note that the values of the electron density in that plane are very small, of the 
order of 10~ 5 e/(a.u.) 3 , so the plane is in the tail region of the electron density distribution. 
Nevertheless the densities clearly reveal the topography of the graphene layer. Electron 
density contours on top of carbon atoms surround other contours representing the large 
hexagonal holes. Densities are larger in the other plane, closer to the plane of nuclei. In 
each plane the density is larger in the positions above carbon atoms and lower above the 
hexagons. A plot that complements this view is given in Figure 5, that gives the electron 
density in a plane perpendicular to the graphene layer through a line containing two adjacent 
carbon atoms, labelled C 1 and C 2 in the figure. Then, points labelled M and X represent 
the midpoint of a carbon-carbon bond and the center of an hexagon respectively. The most 
noticeable feature is the existence of depressions of electron density in the regions above 
the centers of carbon hexagons. These hollow regions are separated by regions of larger 
density that delineate the skeleton of carbon-carbon bonds. In this figure the density of 
the most external contour is p = 1.11 x 10~ 2 e/(a.u.) 3 and the interval between contours 



4 



A/9 = 1.11 x 1(T 2 e/(a.u.) 3 

These observations correlate with the features in Fig. 3, and lead to the following 
interpretation of the potential energy curves. Each curve can be seen as arising from two 
main contributions, one attractive and one repulsive. The attractive contribution is rather 
similar for all the configurations (notice the similarity of the potential energy curves beyond 
6 a.u.) and is mainly due to exchange and correlation effects. Neglecting correlation for the 
purposes of simplicity, the exchange contribution to the total energy is given, in the LDA, 
by the functional 

E L x DA [p\=C x jp{r)U\ (1) 

where C x is a well known negative constant^ In the regime of weakly overlapping densi- 
ties, and assuming no density rearrangements due to the close-shell character of H 2 , the 
contribution of exchange to the interaction energy becomes 

AE X = C X [J [p H2 (r) + p g (r)} I d 3 r - J p H2 (r)l d 3 r - J p g (r)l d 3 r], (2) 

where p g and pu 2 represent the tail densities of the graphene layer and H 2 molecule respec- 
tively. A net "bonding" contribution arises from the nonlinearity of the exchange energy 
functional. On the other hand the sharp repulsive wall is due to the short-range repulsion 
between the close electronic shell of the H 2 molecule and the electron gas of the substrate. 
This contribution is very sensitive to the local electron density sampled by the H 2 molecule 
in its approach to the graphene layer and explains the correlation between the position and 
depth of the different minima in Fig. 3 and the features of the substrate electron density in 
Figs. 4 and 5. Similar arguments explain the physisorption of noble gas atoms on metallic 
surfaces^ and the weak bonding interaction between noble gases £21 At very large separation 
the interaction energy curves should approach the Van der Waals interaction, that is not 
well described, however, by the LDA. 

An interesting point concerns the comparison of the minima of the curves C and D of Fig. 
3. That of curve D is deeper and occurs at a shorter H 2 - graphene separation. The reason 
is that the surfaces of constant electron density of the H 2 molecule have the shape of slightly 
prolate ellipsoids instead of simple spheres. Consequently, for a given distance d between 
the center of mass of H 2 and the graphene plane, the molecule with the perpendicular 
orientation (C) penetrates more deeply into the electronic cloud of the substrate than in 
the parallel orientation (D). In other words, the repulsive wall is reached earlier, that is for 
larger d, in the perpendicular configuration (C). If we consider an electronic density contour 
in H 2 with a value p = 0.018 e/(a.u.) 3 , then the two semiaxes have lengths of 2.07 and 1.71 
a.u. respectively and the difference between these two lengths is 0.36 a.u. This value is in 
qualitative agrement with the difference between the H 2 - graphene separations for the two 
minima of curves C and D, which is 0.20 a.u. This shape effect is usually neglected in the 
phenomenological approaches, that treat H 2 simply as a spherical molecule. 

Figure 6 gives a plot of the charge density difference 

Pdiffij) = Ptot{r) - (p g {r) + pH 2 (r)), (3) 

where Ptot( r ) is the calculated density of the total system, that is the H 2 molecule physisorbed 
in orientation D at a distance of 5 a.u. above the graphene layer, whereas p g + pn 2 is the 
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simple superposition of the densities of the pure graphene layer and H 2 molecule placed also 
in orientation D, 5 a.u. above the graphene layer. That density difference p^jj(r) 
is given in the same plane, perpendicular to the graphene layer, used in Fig. 5. p,ii//( r ) 
has positive and negative regions. The positive region is the area bound by the contour 
labelled P. This region has the shape of two lobes joined by a narrow neck. Contour P 
has a value p^ff = 2.36 x 1CT 5 e/(a.u.) 3 and pdi// increases in this positive region as 
we move towards inner contours in the lobes. The innermost contour shown has a value 
Pdiff — 2.87 x 10~ 4 e/ (a.u.) 3 . The H 2 molecule sits above the neck, so the figure reveals that 
the repulsive interaction produced by the close electronic shell of H 2 pushes some charge 
from the region immediately below the molecule (the neck region) to form the lobes of 
positive Pdiff (r)- This displacement of electronic charge is nevertheless quantitatively very 
small. Notice that p g takes values between 1.6 x and 4.1 x 10~ 3 e/(a.u.) 3 in a plane 
3 a.u. above the graphene layer, while p^iff has values of the order 1CT 5 — 1CT 4 e/(a.u.) 3 
in the same plane. The smallness of pdiff justifies the argument given in eq. (2) for the 
attractive exchange-correlation contribution to the interaction potential. 

The static calculations discussed above have been complemented with dynamical simula- 
tions in which the H 2 molecule was initially placed in different orientations at distances of 4 - 
6 a.u. from the graphene layer and was left to evolve under the influence of the forces on the 
H atoms. The H 2 bondlength was allowed to adjust in the process. The simulations confirm 
the results of the static calculations, in the sense that the H 2 molecules end up in positions 
above the center of an hexagon at the end of the simulations. The binding energies and 
H 2 - graphene layer distances practically coincide with those in Table II. Marginally small 
differences in separation or binding energy are due to very small changes of the bondlength 
of H 2 , always smaller than 0.3%. The result of one of the simulations is worth to be men- 
tioned. A configuration intermediate between those labelled C and D above was obtained: 
the center of mass of the molecule was 5.10 a.u. above the center of a carbon hexagon, with 
the molecular axis forming an angle of about 30° with the graphene plane. The binding 
energy in this new configuration was only 1 meV larger than in the parallel configuration D. 

In summary, the picture arising from the calculations is rather clear. The H 2 molecules 
prefer the hollow sites above the centers of carbon hexagons where the background electron 
density is lower than in channels on top of the skeleton of carbon-carbon bonds. The 
exchange-correlation contribution provides the weak attraction responsible for physisorption, 
but the preferred distance of approach is determined by the repulsive part of the interaction 
potential. That repulsive contribution is due to the close-shell electronic structure of H 2 . We 
have performed static calculations of the barrier for the diffusion of a molecule, initially in 
the parallel configuration D at the preferred distance of 5.07 a.u. above the graphene plane, 
to an equivalent configuration D above an adjacent hexagon. The initial configuration of 
the molecule, with its axis perpendicular to two parallel carbon-carbon bonds, can be seen 
in the bottom panel of Fig. 1. The molecule was then forced to follow a path across one 
of those bonds, allowing for the reorientation of the molecular axis at each step in order 
to minimize the energy of the system. Although the molecule begins with the axis parallel 
to the graphene plane, the orientation of the axis changes as the molecule approaches the 
carbon-carbon bond. In fact, when the center of mass of the molecule is precisely above that 
bond, the molecular axis becomes perpendicular to the graphene plane, that is the molecule 
adopts configuration B, as indicated also in Fig. I. The energy difference between this 
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saddle configuration and the starting one gives a calculated diffusional barrier of 14 meV. 
A temperature of 163 K is enough to surpass this barrier. 

The conclusions from the calculations are, in our view, general enough that one can 
make some extrapolations to the case of adsorption of H 2 by carbon nanotubes. When 
adsorption occurs on the outside wall of an isolated nanotube, the predictions of Fig. 3 will 
be valid, with a minor influence of the nanotube curvature. If the tubes form a parallel 
bundle and we consider the interstitial channels between tubes, the effects seen in Figure 3 
will be smoothed out because of the addition of contributions of different graphitic surfaces 
non in registry. Addition of these contributions will give rise to an interstitial channel with 
a potential energy nearly independent of z, if we call z the direction parallel to the tube 
axis. Finally, the same smoothing effect will occur in the inner channel of a tube if the 
tube diameter is not large. In summary we predict very easy diffusion of the H 2 molecule 
in arrangements of parallel tubes along the direction parallel to the tube axis, both inside 
the tube cavity and in the interstitial channels. 

The present adsorption results can be partialy compared with those of Stan and ColeS 
They considered the H 2 molecule as a spherically symmetric entity and calculated the ad- 
sorption potential inside zigzag (13,0) nanotubes (radius = 9.62 a.u.) based on a sum of 
isotropic Lennard- Jones interactions between the molecule and the carbon atoms of the 
tube. Our calculation and that of Stan and Cole agree in that the smallest binding energy 
is obtained for the H 2 upon one carbon atom and the largest one for the H 2 upon the center 
of the hexagon of carbon atoms. However Stan and Cole do not distinguish between parallel 
and perpendicular orientations because they considered an spherical molecule. Their Fig. 1 
shows a binding energy about 0.079 eV for adsorption in front of the center of an hexagon 
of carbon atoms and that the equilibrium distance between the molecule and the nanotube 
wall is 5.7 a.u. This distance is consistent but a little larger than those reported in our 
Table II. On the other hand, the value 0.079 eV for the binding energy is also consistent 
with the binding energies in Table II. Notice, however, that the binding energy for a tube of 
larger radius, or for a planar graphene sheet, will be a little smaller because the curvature 
of the tube increases the number of nearest neighbor carbon atoms. In fact, Wang and 
Johnson % calculated an adsorption binding energy near 0.050 eV for molecular hydrogen in 
an idealized carbon slit pore with a pore width of 17.4 a.u. 



IV. CONCLUSIONS 

By performing DFT calculations we confirm that physisorption of H 2 on graphitic layers 
is possible. The differences between the binding energies corresponding to different positions 
(on top of carbon atoms, on top of carbon-carbon bonds, on top of hexagonal holes) are small, 
and the diffusional barriers are also small, so easy diffusion is expected at low temperature. 
The nonsphericity of the H 2 molecule has some influence on the preferred orientation of the 
molecular axis with respect to the graphene plane. These small effects associated to different 
positions and orientations of the physisorbed molecule are expected to average out inside 
carbon nanotubes or in the interstitial channels in parallel arrays of carbon nanotubes. 
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TABLES 

TABLE I. Calculated energy (Ry) of the graphene layer for several layer-layer distances. The 
energies are calculated for a plane wave cut-off energy of 40 Ry. 



Layer-layer distance (a.u.) 


Energy (per atom) 


20 


-11.4234 


25 


-11.4235 


30 


-11.4235 


35 


-11.4235 



TABLE II. Binding energy (eV) and equilibrium distance (a.u.) for H2 physisorbed on a 
graphene layer. A, B and C correspond to the configurations in which the molecular axis is 
perpendicular to the graphene plane and the molecule is on top of: a carbon atom (A), the midpoint 
of a carbon-carbon bond (B), the center of an hexagon (C). In configuration D the molecule is on 
top of the center of an hexagon with the molecular axis parallel to the graphene plane. 





A 


B 


C 


D 


Binding energy 


0.070 


0.072 


0.083 


0.086 


Distance 


5.50 


5.49 


5.25 


5.07 
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FIGURES 



FIG. 1. Top panel gives a fragment of the graphene layer showing the eight carbon atoms in 
the unit cell, represented by large spheres. Bottom panel shows the three adsorption configurations 
with the molecular axis perpendicular to the graphene plane. These have the H2 molecule above 
a carbon atom (A), above the midpoint of a carbon-carbon bond (B) and above the center of an 
hexagon (C). Also shown is an adsorption configuration (D) with the molecular axis parallel to the 
graphene plane and the molecule above the center of an hexagon. 

FIG. 2. Comparison of potential energy curves for the parallel approach of H2 to the graphene 
layer upon the center of an hexagon of carbon atoms. The curves were obtained using supercells 
such that the graphene layers are separated by 30 a.u. (circles) or 20 a.u. (crosses). 

FIG. 3. Potential energy curves for the aproach of H2 to the graphene layer in four different 
configurations. The axis of the molecule is perpendicular (A, B, C) or parallel (D) to the graphene 
layer. In the former orientation the molecule is above a carbon atom (A), above the center of a 
C-C bond (B), and above the center of an hexagon (C). In the parallel orientation (D) the molecule 
is above the center of an hexagon. 

FIG. 4. (a) Contours of constant electron density p of a pure graphene layer in a plane 5 a.u. 
above the plane of the carbon nuclei, p = 8.64 x 10 -5 e/(a.u.) 3 in the innermost contours above 
carbon atoms, p = 6.08 x 10~ 5 e/ (a.u.) 3 in the innermost contours above the large hexagonal holes. 
Densities decrease monotonously between those two contours with an interval Ap = 0.18 x 10~ 5 
e/(a.u.) 3 . (b) Contours in a plane 3 a.u. above the plane of carbon nuclei, p = 3.92 x 10~ 3 e/ (a.u.) 3 
in the innermost contours above carbon atoms, p = 1.72 x 10 -3 e/(a.u.) 3 in the innermost contours 
above the large hexagonal holes. Densities decrease monotonously between those two contours with 
an interval Ap = 0.16 x 10~ 3 e/(a.u.) 3 

FIG. 5. Contours of constant electron density of pure graphene in a plane perpendicular to the 
graphene layer, going through a line containing two adjacent carbon atoms, labelled C\ and C2. 
Symbols M and X indicate the mid-point of a carbon-carbon bond and the center of an hexagon, 
respectively. The outermost contour plotted is p = 1.11 x 10 -2 e/(a.u.) 3 and the interval between 
contours Ap = 1.11 x 10~ 2 e/(a.u.) 3 

FIG. 6. Charge density difference Pdiff = Ptot ~ Pg ~ Ph 2 f° r ^2 physisorbed 5 a.u. above the 
graphene layer. The plane of the plot and the symbols CI, C2, M and X are the same as in Fig. 
5. Contour labelled P has a value 2.36 x 10~ 5 e/(a.u.) 3 and encloses the region of positive Pdiff- 
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